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Abstract 

We  discuss  how  to  describe  the  Markov  chain  underlying  a  generalized  stochastic  Petri  net  using  Kro¬ 
necker  operators  on  smaller  matrices.  We  extend  previous  approaches  by  allowing  both  an  extensive 
type  of  marking-dependent  behavior  for  the  transitions  and  the  presence  of  immediate  synchronizations. 
The  derivation  of  the  results  is  thoroughly  formalized,  including  the  use  of  Kronecker  operators  in  the 
treatment  of  the  vanishing  markings  and  the  computation  of  impulse-based  reward  measures.  We  use 
our  techniques  to  analyze  a  model  whose  solution  using  conventional  methods  would  fail  because  of  the 
state-space  explosion.  In  the  conclusion,  we  point  out  ideas  to  parallelize  our  approach. 
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1  Introduction 


Generalized  stochastic  Petri  nets  (GSPNs)  [2,  3,  6]  are  ideally  suited  to  model  a  large  class  of  performance 
and  reliability  problems,  but  their  numerical  analysis  requires  the  solution  of  a  very  large  continuous¬ 
time  Markov  chain  (CTMC).  The  size  of  the  transition  rate  matrix  R  for  the  CTMC  is  the  main  obstacle, 
since  its  memory  requirements  can  easily  exceed  the  capacity  of  today’s  (and  tomorrow’s)  machines  even 
when  sparse  storage  techniques  are  employed. 

A  possible  approach  to  this  problem  is  to  store  R  implicitly.  Plateau  [14]  proposed  the  use  of 
Kronecker  operators  for  the  description  of  the  transition  rate  matrix  of  a  structured  model  composed 
of  a  set  of  “synchronized”  stochastic  automata,  a  subclass  of  GSPNs.  Buchholz  [4]  used  a  similar  idea 
for  Markovian  closed  asynchronous  queueing  networks,  and  Takahashi  [17]  used  it  for  open  queueing 
networks  with  communication  blocking.  Donatelli  [10,  11]  adapted  the  approach  to  GSPNs,  defining 
first  the  “superposed  stochastic  automata”,  then  the  “superposed  GSPNs”.  Later,  Buchholz  [5]  applied 
the  concept  to  the  special  case  of  marked  graphs  and  Kemper  [12]  addressed  some  of  the  problems  in 
[11],  extending  the  applicability  of  the  results. 

These  approaches  have  in  common  a  decomposition  of  the  model  into  a  set  of  submodels,  so  that 
the  state  space  of  the  CTMC  underlying  the  entire  model  is  a  subset  of  the  cross-product  of  the  state 
spaces  of  the  CTMCs  underlying  each  submodel.  This  implies  that  the  transition  rate  for  the  entire 
model  can  be  described  using  Kronecker  operators  on  smaller  matrices. 

Focusing  on  GSPNs,  the  result  of  decomposing  a  GSPN  is  a  set  of  largely  independent  sub-GSPNs, 
but  some  transitions  will  be  shared  by  multiple  sub-GSPNs,  to  model  interactions  among  them.  If  a 
transition  tj  is  shared  by  two  sub-GSPNs  A\  and  A2,  and  tj  has  input  and  output  places  in  both  of 
them,  this  models  a  synchronization.  Aj  and  A2  “wait  for  each  other”  and  the  event  corresponding  to 
tj  occurs  only  when  they  are  both  “ready”.  An  alternative  case  arises  when  tj  has  its  input  places  in  A] 
and  its  output  places  in  A2.  This  describes  an  asynchronous  (and  asymmetric)  communication,  since  A2 
must  “wait  for  permission”  from  A,.  However,  if  an  output  place  of  tj  in  A2  has  a  capacity  defined  for 
it,  Aj  will  wait  for  A2  as  well.  One  of  the  contributions  of  our  work  is  to  present  a  unified  framework  for 
all  these  types  of  interactions.  Indeed,  we  do  not  assume  that  the  net  possesses  any  particular  structure. 

The  solution  of  a  decomposed  GSPN  is  then  based  on  the  following  idea  [11],  Let  n,  be  the  number 
of  states  in  Sf,  the  state  space  for  the  CTMC  underlying  the  t-th  sub-GSPN,  i  =  S'T  = 

S j.  x  . . .  x  Sj  3  St  is  the  “potential  state  space”  for  the  model,  usually  (much)  larger  than  the  actual 
state  space  St,  so  that  a  transition  rate  matrix  R'  is  defined  using  Kronecker  algebra: 

R/  =  0R’+  £  (g)R*-*, 

«=1  (jgT*  »=1 

where  R’  describes  the  local  transitions  for  the  i-tli  sub-GSPN  (including  the  effect  of  immediate 
transitions),  T*  is  the  set  of  synchronizing  transitions,  which  must  be  timed,  and  the  “corrective  matrix” 
R*’-7  describes  the  effect  of  tj  on  the  *-th  sub-GSPN. 

4  lie  actual  tiansition  rate  matrix  R  can  be  obtained  from  R/  by  eliminating  the  rows  and  columns 
corresponding  to  the  unreachable  states  in  Sf  \  ST,  but  this  requires  an  additional  overhead,  since 
the  composition  of  a  marking  does  not  directly  indicate  whether  it  is  reachable  or  not.  Hence,  an 
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alternative  approach  was  initially  suggested.  The  power  or  Jacobi  method  is  used  to  compute  the 
steady-state  solution  in  a  vector  ir'  of  size  |<Sj|.  By  assigning  a  nonzero  initial  probability  only  to 
markings  in  St,  the  solution  7r'  is  guaranteed  to  be  zero  for  any  marking  m  €  S'T  \  St-  This  simplifies 
the  algorithm,  since  m  can  now  be  interpreted  as  the  mixed-base  integer  index  of  the  corresponding 
entries  in  R'  and  7r',  but  the  memory  requirement  might  be  excessive  when  |<S^|  |<Sr|. 

To  reduce  the  impact  of  unreachable  markings,  Kemper  [12]  proposed  a  technique  that  only  requires 
a  probability  vector  7r  of  size  |<Sr|  •  In  the  numerical  iterations,  for  each  m  6  St,  each  entry  R[uu  >  0 
is  generated  (this  implies  n  6  St)  and,  given  n,  its  index  k  in  7r,  or  its  lexicographic  position  in  St,  is 
computed  in  0(log  |Sr|)  operations,  using  a  binary  search. 

We  unify  previous  work  by  offering  a  thorough  discussion  of  the  structure  of  the  underlying  CTMC, 
including  the  management  of  immediate  transitions  and  vanishing  markings.  Our  formalism  is  more 
general  than  those  assumed  in  [4,  5,  10,  11,  12],  since  we  allow  for  marking-dependent  arc  cardinali¬ 
ties  and  rates,  subject  to  certain  restrictions,  hence  our  results  include  those  previously  mentioned  as 
special  cases.  Then,  using  an  approach  based  on  discrete-time  Markov  chains  (DTMCs),  we  also  re¬ 
move  the  main  restriction  previously  imposed  on  the  decomposition  of  the  GSPN:  we  allow  immediate 
synchronizing  transitions.  Finally,  we  consider  a  reward  structure  defined  on  the  GSPN,  and  we  show 
how  to  compute  the  expected  reward  in  steady-state  in  the  Kronecker  framework.  This  is  of  particular 
importance  for  impulse  rewards  associated  with  immediate  transitions,  whose  firing  are  only  implicitly 
represented  in  R. 

The  paper  is  structured  as  follows.  Section  2  describes  the  notation  used  and  recalls  the  main 
concepts  of  Kronecker  algebra,  GSPNs,  Markov  chains,  and  rewards.  Section  3  presents  the  expression 
for  the  transition  rate  matrix  of  the  CTMC  underlying  a  generic  GSPN,  provided  its  transitions  satisfy 
certain  restrictions  on  the  type  of  marking-dependency.  The  result  is  quite  general,  but  not  directly 
applicable,  since  it  requires  one  to  compute  the  inverse  of  a  matrix  described  as  the  sum  of  Kronecker 
products.  However,  Sections  4  and  5  use  it  to  derive  computationally  effective  expressions  for  GSPNs 
with  timed  and  immediate  synchronizing  transitions,  respectively.  Implementation  and  application  of 
these  results  are  shown  in  Section  6,  including  detailed  information  about  computation  time  and  memory 
requirements.  Section  7  contains  a  summary  and  discusses  future  extensions,  including  distributed 
implementations  and  approximate  solutions. 


2  Notation  and  definitions 

Except  for  W,  the  sets  of  natural  numbers,  {0,1,2...},  and  2R,  the  set  of  real  numbers,  all  sets  are 
denoted  by  upper  case  calligraphic  letters  (e.g.,  A );  vectors  and  matrices  are  denoted  by  lower  and 
upper  case  bold  letters,  respectively  (e.g.,  a,  A);  their  entries  are  denoted  by  subscripts  (e.g.,  a*,  A 
a  set  of  indices  can  be  used  instead  of  a  single  index,  for  example,  Axy  denotes  the  submatrix  of 
A  corresponding  to  set  of  rows  X  and  the  set  of  columns  y .  Superscripts  denote  families  of  related 
quantities  (e.g.,  A1,  A2).  0XXy  and  lxXy  denote  matrices  with  x  rows  and  y  columns,  having  all  entries 
equal  to  0  or  1,  respectively,  while  Ix  denotes  the  identity  matrix  of  size  x  x  x;  the  dimensions  of  these 
matrices  are  omitted  when  they  are  clear  from  the  context.  Given  a  vector  a,  diag(a)  is  a  square  matrix 
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having  a  on  the  diagonal  and  zeros  elsewhere.  Given  annxn  matrix  A,  rowsum(A)  =  diag(A  ■  lnxl) 
is  a  matrix  having  the  diagonal  equal  to  the  sums  of  the  entries  on  each  row  of  A,  and  zeros  elsewhere, 
while  6(A)  is  a  matrix  having  the  same  nonzero  pattern  as  A,  but  with  entries  equal  to  either  0  or  1. 

2.1  Kronecker  algebra 

We  recall  the  definition  of  the  Kronecker  product  B  =  Ak  of  K  matrices  Ak  €  jRn*xm*,  using  the 
convention  that  the  rows  and  columns  of  both  B  and  the  matrices  Ak  are  indexed  starting  at  0.  The 
generic  element  of  B  G  fRrijt=i  n*xIlLi is 

B(...((»l)n2+«2)n3-)njt+»fcl(...((j1)m2+i2)m3-)»»*+it  =  ^:i  ,J,  A?2)J2  *  *  *  ,JK 

with  0  <  ik  <  nk  and  0  <  jk  <  m*,  for  k  =  1 Assuming  a  mixed-base  numbering  scheme 
so  that  the  tuple  (/,,  /2, . . .  Ik)  corresponds  to  row  (,..((/j)n2  +  /2)n3  •  •  -)nk  +  lk  or  column  (...((/i)j7i2  + 
h)mz  •  •  respectively,  we  will  also  write  the  above  quantity,  more  succinctly,  as  B 

The  Kronecker  sum  ©*L,  Ak  of  I<  square  matrices  A*  €  RnkXllk  is  defined  as 

K  K 

0  Ak  =  £  I„,  ®  •  •  •  <g>  Injt_,  ®  Ak  ®  Injb+1  •  •  •  ®  I„A.. 

*=1  Jk=l 

2.2  Generalized  stochastic  Petri  nets 

A  generalized  stochastic  Petri  net  (GSPN)  is  a  tuple  (V.T,!,  C~,C+,m°,w),  where: 

•  P  =  {Pi» •  •  •  ,P\v\)  is  ^  finite  set  of  places ,  drawn  as  circles  in  the  graphical  representation  of  the 
GSPN.  A  non-negative  integer  vector  m  6  called  marking  describes  the  number  of  tokens  in 
each  place.  Given  a  place  p,  €  V,  m,  is  the  number  of  tokens  in  pi  for  marking  m. 

•  T  =  {ti , . . . ,  ^|Tj}  is  a  finite  set  of  transitions ,  V  fl  T  =  0. 

•  3T  C  T  is  the  subset  of  immediate  transitions,  drawn  as  thin  bars,  while  X  =  T  \  J  are  the  timed 
transitions ,  drawn  as  rectangles.  The  firing  time  of  immediate  transitions  is  zero,  while  that  of 
timed  transitions  is  exponentially  distributed. 

•  C“  and  C+  are  incidence  matrices  of  size  \P\  X  |T|.  Their  elements  are  functions  from  IN ^  to 
IN-  Qj(m)  and  C+^m)  denote  the  marking-dependent  integer  cardinality  assigned  to  the  input 
arc  from  7^  to  tj  and  the  output  arc  from  tj  to  pi  respectively.  In  the  graph,  these  arcs  are  drawn 
using  an  arrowhead  pointing  to  the  destination  if  their  cardinality  is  not  identically  equal  to  zero. 
The  cardinality  function  is  indicated  on  the  arc  unless  it  is  identically  equal  to  one. 

•  m°  is  the  initial  marking .  In  the  graph,  the  value  of  m-  is  written  inside  place  pt,  if  positive. 

•  For  any  tj  €  T,  Wj  is  a  function  from  to  1R.  Wj(m)  is  the  weight  associated  with  transition 

tj  in  marking  m.  According  to  whether  tj  is  immediate  or  timed,  this  weight  represents  an 
(unnormalized)  firing  probability ,  or  a  firing  rate. 
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A  transition  tj  6  T  has  concession  in  marking  m  iff 

Vp,-  €  V,  C";(m)  <  m„  or  ^(m)  <  m. 

If  any  immediate  transition  has  concession  in  m,  it  is  enabled,  and  m  is  said  to  be  vanishing.  Otherwise, 
m  is  said  to  be  tangible  and  any  timed  transition  tj  with  concession  is  enabled  in  m.  In  other  words,  a 
timed  transition  is  not  enabled  in  a  vanishing  marking  even  if  it  has  concession. 

Some  definitions  of  GSPNs  allow  one  to  disable  a  transition  tj  with  concession  in  m  by  specifying  a 
zero  weight  Wj(m)  for  it,  or  by  introducing  inhibitor  arcs,  drawn  with  a  circle  instead  of  an  arrowhead. 
In  a  marking  m,  an  inhibitor  arc  from  place  p,-  to  transition  tj  with  cardinality  c(m)  disables  tj  if  m,  > 
c(m).  Since  these  behaviors  can  be  represented  by  an  appropriate  definition  of  input  arc  cardinalities 
in  our  formalism,  we  assume,  without  loss  of  generality,  that  Wj(m)  >  0  if  tj  is  enabled  in  m,  and  we 
use  inhibitor  arcs  in  our  models  merely  as  a  shorthand. 

Let  £(m)  denote  the  set  of  enabled  transitions  in  marking  m.  An  enabled  transition  tj  firing  in 
marking  m  yields  a  new  marking  n  such  that 

Vp,  €  V,  n,  =  m,  -  C“,(m)  +  C+ (m)  =  m,  +  CXii(m)  (or  n  =  m  +  C7.,i(m)), 

where  C  =  C+  —  C"  is  the  incidence  matrix  of  the  GSPN.  We  can  also  write  m-^n  to  express  that  tj 
has  concession  in  m  and  that  n  is  obtained  from  m  by  firing  tj,  regardless  of  whether  tj  6  £(m)  or  not 
(tj  is  not  enabled  if  it  is  timed  and  m  is  vanishing,  or  if  w3(m)  =  0). 

The  firing  probability  of  a  transition  tj  enabled  in  marking  m  is 

wj(m)  (j) 

If  m  is  tangible,  this  corresponds  to  a  race  between  the  exponentially  distributed  firing  times  of  the 
enabled  transitions,  with  rates  given  by  w.  In  a  vanishing  marking,  instead,  weights  define  a  probabilistic 
choice,  since  the  race  model  does  not  specify  how  to  choose  which  transition  to  fire  next  when  multiple 
enabled  transitions  have  the  same  zero  firing  time. 

2.3  Reachability  set 

The  reachability  set  S  is  defined  as  the  set  of  markings  reachable  from  the  initial  marking  m°  by  firing 
any  sequence  of  enabled  transitions.  Formally,  <S  is  the  smallest  subset  of  containing  m°  and  such 
that  m  €  S,  tj  €  £(m),  and  m-^n  imply  n  £  S.  Fig.  1  shows  the  skeleton  of  an  algorithm  to  build  the 
set  of  reachable  markings  S  (which  we  assume  finite  from  now  on).  Particular  care  must  be  placed  on 
the  implementation  of  statement  9,  since  the  size  of  the  set  S  to  be  searched  is  very  large  in  practice. 
Efficient  methods  include  hashing  or  balanced  search  trees  (e.g.,  AVL  trees  [18]).  While  not  explicitly 
stated  in  the  algorithm,  S  should  be  stored  as  the  union  of  two  disjoint  sets,  St  and  <Sk,  corresponding 
to  the  tangible  and  vanishing  markings,  respectively. 

A  function  'P  assigns  an  index  to  each  reachable  marking,  according  to  a  lexicographic  order,  indi¬ 
cated  by 

:  S  — >  {0, . . . ,  |<S|  —  1}  such  that  'P(m)  >  ^(n)  ■$=$>  m  >-  n. 
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Algorithm  BuildRS  (input:  (V>T,I,C  ,C+,m°,w);  output:  <S); 

1.  S  <—  {m0};  /*  S  contains  the  markings  found  so  far  */ 

2.  U  <—  {m0};  /*U  C  S  contains  the  found  but  unexplored  markings  */ 

3.  while  U  ^  0  do 

4.  “choose  a  marking  m  from  ZY”; 

5.  U  <r-U  \  {m}; 

6.  “compute  £(m)”; 

7.  for  each  j  6  £(ni)  do 

8.  n  <—  m  +  Cpj(m); 

9.  if  n  £  S  then 

10.  IA  < —  U  U  {n}; 

11.  S  <—  S  U  {n}; 

12.  end  if 

13.  end  for 

14.  end  while 


Figure  1:  Algorithm  BuildRS 

If  an  AVL  tree  is  used,  lP  can  be  precomputed  with  a  simple  preorder  visit  of  the  tree,  and  its  value  can 
be  stored  in  the  nodes  of  the  tree.  Then,  given  m  £  the  value  k  =  ty(m)  can  be  found  in  0(log  |<S|) 
operations  using  the  AVL  tree  augmented  with  this  additional  information  ('l'(m)  =  “undefined”  for 
any  m  ^  5).  The  restrictions  of  to  the  tangible  and  vanishing  markings,  can  be  defined  accordingly: 
:  St  — ►  {0, . . . ,  }St\  —  1}  and  tyy  :  Sy  — »  {0, . . . ,  |<SV|  —  1}- 

In  the  following,  with  a  slight  overloading  in  the  notation,  we  use  a  marking  m  to  index  data 
structures  (vectors  and  matrices)  referring  to  <S,  Sr ,  or  Sy.  Strictly  speaking,  we  should  use  instead 
'l'(m),  m),  and  'l'v'(m),  respectively,  but  this  would  make  the  expressions  excessively  cumbersome. 

Nevertheless,  it  is  important  to  stress  this  fundamental  difference  from  a  computational  point  of  view; 
finding  the  index  of  a  marking  is  a  potentiaLsource  of  additional  complexity  in  any  structured  approach. 

2.4  Underlying  continuous-time  Markov  chain  and  rewards 

We  focus  on  the  steady-state  analysis  of  the  continuous-time  Markov  chain  (CTMC)  underlying  a  GSPN, 
described  by  the  infinitesimal  generator  matrix  Q  €  which  we  assume  ergodic  (and  finite, 

since  S  is  finite): 

Q  =  R  —  A  —  R t,t  +  Rr,v(I|5v|  —  U v7v)^Vyj  -  A,  (2) 

where  R  is  the  transition  rate  matrix,  A  =  rowsum(R),  and  Rr.r  and  R ry  (U yj  and  U^v)  describe 
the  rates  (probabilities)  of  going  from  tangible  (vanishing)  markings  to  tangible  or  vanishing  markings, 
respectively.  The  entry  of  Rt,t  (Rr,v)  corresponding  to  the  row  for  m  and  the  column  for  n  describes 
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the  rate  from  m  E  Sr  to  n  E  Sr  (n  E  Sy): 


Y  wj(m). 

tj 

tj€£(m),m — rii 

The  entries  of  Uy,v  and  U v,t  are  defined  analogously,  using  the  firing  probability  of  immediate  transition 
tj,  given  by  (1),  instead  of  the  weight  Wj(m).  For  a  discussion  of  how  to  generate  Q  in  practice,  see 
[2,6]. 

We  observe  that  /lm>m  =  X)tjg£(m)  wj(m)  is  then  the  total  rate  leaving  marking  m  €  St,  and  it 
equals  the  inverse  of  hm,  the  expected  holding  time  in  m. 

Let  7rm  be  the  steady-state  probability  of  a  tangible  marking  m  (vanishing  markings  have  zero 
probability).  Then,  the  steady-state  probability  (row)  vector  7r  £  IR)St\  satisfies  the  balance  equation 


7r  •  Q  =  01x|5t|  subject  to  the  normalization  ir  •  l|sr|Xi  =  1. 


(3) 


We  can  specify  a  quantity  of  interest  for  the  GSPN  using  a  reward  structure  (/>,  r),  where  />(m)  is 
the  reward  rate  gained  while  the  GSPN  is  in  marking  m,  and  rj(m)  is  the  reward  impulse  gained  when 
transition  tj  E  T  fires  in  marking  m.  The  expected  reward  rate  in  steady  state  is  then 

Y  TTmp(m)  +  J2  Y  $i,mri( m),  (4) 

m eST  m eS  qef(m) 


where  4>Jim  is  the  rate  at  which  transition  tj  fires  in  steady  state  in  marking  m.  If  we  let  (f)  £  IR)S\  be 
the  vector  describing  the  rate  at  which  each  marking  is  entered  in  steady  state,  is  obtained  as 


j 


wi(m) 


111 


.q6£(.n)W/(m) 

For  m  £  St,  4>m  =  w/(m)  =  7rmAmm.  For  m  €  Sy,  instead, 

4* m  =  Yj  7r‘1  ‘ 

ii£St 


where,  for  n  £  St  and  m  £  Sy,  the  corresponding  entry  of  matrix 


F  =  Rt,v(I|5h  -  Uv.v)"1  €  JR15*'1*15"1  (5) 

describes  the  rate  at  which  a  vanishing  marking  m  is  entered  after  leaving  a  tangible  marking  n  and 
before  reaching  the  next  tangible  marking.  If  no  reward  impulse  is  defined  for  immediate  transitions, 
then  Eq.  (4)  reduces  to 

Y  ( Mm)  +  Y  wi(m)rj(m) 

m  esT  \  t,ee(m) 

In  this  work,  wc;  consider  the  structure  of  both  Q  and  F,  and  the  computation  of  7T  and  cf). 
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3  Kronecker  expression  for  the  CTMC  underlying  a  GSPN 

We  now  show  how  the  ideas  in  [4,  5,  10,  11,  12]  can  be  applied  to  individual  places,  not  just  to  sub- 
GSPNs.  Our  goal  is  to  clarify  the  relationship  between  Kronecker  algebra  and  GSPNs,  while  relaxing 
several  important  restrictions  on  the  type  of  interactions.  Later  we  will  merge  individual  places  into 
“macroplaces”,  which  corresponds  to  the  notion  of  sub-GSPNs. 


3.1  Using  the  state  spaces  of  individual  places 

The  first  extension  regards  the  type  of  marking-dependency  allowed  in  the  GSPN.  We  allow  the  weight 
of  a  transition  to  be  expressed  as  the  product  of  “local  effects”  due  to  the  number  of  tokens  in  each 
place: 

Vm  €  <S,  Vi,  E  T,  C^..  <  m  =>  w^m)  =  iu*  •  (6) 

Piev 

where  Wj  can  be  interpreted  as  a  constant  “reference”  weight,  while  the  values  Wij  are  dimensionless 
scaling  functions  [1].  This  “independence  of  effects”  in  the  marking  dependence  implies,  for  example, 
that  if  markings  m  and  n  differ  only  in  the  number  of  tokens  in  and  if  ij  is  enabled  in  both, 
Wj(n)  =  Wj(m)  •  Wij(ni)/witj(mi).  If  a  weight  w j  does  not  depend  on  the  number  of  tokens  in  p,-,  we 
assume,  without  loss  of  generality,  that  Wij  is  identically  equal  to  one.  Note  that  we  do  not  require  the 
weight  Wj(m)  of  a  timed  transition  tj  in  a  vanishing  marking  m  to  be  zero.  Doing  so  would  make  the 
specification  of  w  for  a  given  GSPN  more  difficult  in  practice,  and  is  not  required  by  our  approach. 
Analogously,  the  dependence  of  the  matrices  C”  and  C+,  hence  C,  is  assumed  to  be  of  the  form 

Vm  €  5,  €  Vyt3  E  T,  C^.(rn)  =  (7) 

where  is  one  of  or  nothing,  and  fifj  is  a  function  from  IN  to  IN. 

We  can  now  state  a  theorem  expressing  the  matrices  Q  and  F  of  the  CTMC  underlying  a  GSPN  in 
terms  of  smaller  matrices  related  to  each  place-transition  pair. 


Theorem  3.1  Consider  a  GSPN  with  finite  reachability  set  S  satisfying  Eq.  (6)  and  (7), 
and  let  n,—  1  be  the  bound  of  place  pi  E  V,  that  is,  for  any  m  E  <S,  m*  E  {0, 1, . .  .n,  —  1}  =  S\ 
Define 

r'=  £  .  0  w'’j'  w,j,  (8) 

tj€X  v,£7>  tj£l  p,€V 

where  W,J  is  a  square  matrix  of  size  n,  x  whose  entry  in  position  (r,c),  for  r,c  €  S\  is 
given  by 


W*'J'(r,  c) 


wi,j{r)  if  r  >  Pitj(r )  and  c  =  r  +  (3itj(r) 
0  otherwise 


Also,  define 


A'  =  rowsum(R')  —  wj  '  ^  rowsum(W*,•,)  =  ^  Wj  •  (^)  (9) 

tj€.v  P,ev  tjex  p,€T 
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r  =  rowsum(U')  =  £  w*  •  (g)  rowsum(Wi*>)  =  £  u$  •  (g)  riJ  (10) 

iyel  p;€7>  qeX  P.-6T> 

and  T'  =  I|5-|  —  S(r').  Then,  the  matrix 

Q'  =  R'  (l-(T'  +  r)"'  •U')-'  -  A'  (11) 

satisfies  Q  =  Q^,  St  and  F  =  Q^,^,  where  Q  and  F  have  the  meaning  defined  in  Eq.  (2) 
and  (5). 

Proof:  Matrices  with  superscript  have  row  and  column  set  S1  =  {o,i,...(iw»0  - !}’ 
or  S1  x  •  •  •  x  if  we  identify  a  tuple  with  its  mixed-base  value.  In  the  following,  however, 
we  partition  matrices  and  permute  their  rows  and  columns  so  that  the  markings  appear  in 
lexicographic  order  within  the  sets  St,  Sy,  and  S'  \  S.  This  is  for  illustration  purposes  only. 

First,  we  prove  that 


Rm.n  =  ( E  w"i  ■  ®  W'J)  =  E  ■  II  =  E  w>(m)'  (12) 

V,5.v 

Let’s  consider  the  contribution  to  this  value  for  each  timed  tj  £  X,  by  doing  a  case  analysis: 

1.  If  m-^n,  the  contribution  should  be  Wj(m).  Indeed,  for  all  p,  €  V,  =  u>tij(m,), 

hence  the  contribution  of  tj  is 

wj  ■  II  wiAmi)  =  Wj(m). 

p>ev 


2.  If  tj  does  not  have  concession  in  m  the  contribution  should  be  zero.  Indeed,  there  must 
exist  a  place  p,  such  that  m,  <  C£,-(m)  =  fi~j(mi).  This  implies  WJ^.  n.  =  0,  and  the 
contribution  of  tj  is  w*  •  rww£,ni  =  o. 

3.  If  m^n'  7^  n,  the  contribution  of  tj  should  be  zero  as  well.  Indeed,  there  must  exist 

a  place  p,  such  that  n,  ^  m,  -f  C,  j(m)  =  m,-  -f  lienee,  W^.  n.  =  0,  and 

w*j  ■  npi67>Win„u.  =  0- 

Thus,  the  contribution  of  each  transition  in  the  summation  is  correct.  An  analogous  argu¬ 
ment  allows  us  to  show  that 

um,n  =  Wj(m).  (13) 

ty€  J.m-i-n 


From  Eq.  (12)  and  (13),  we  can  conclude  that  =  0,  U^,  =  0,  R^r  <s,^5  =  0,  and 

that  the  matrices  R t,t,  Rr.v,  Uv.r,  and  Uv.v  for  the  underlying  GSPN,  with  their  rows 
and  columns  ordered  according  to  VI/,  can  be  expressed  as: 


Rt,t  —  RsTlsr  Rr.v  -  R'st,sv  ^v,t  —  r'Sv>Sv  ■  U'Sv  St  Uv'.v  —  r'sv,sv  ' 


Sy  ,S\/ 

(14) 
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(the  normalization  r'2ltsv  IS  required  because  the  weights  of  the  immediate  transitions 
enabled  in  a  vanishing  marking  are  not  required  to  sum  to  one,  while  the  entries  in  U y? 
and  U yy  are  probabilities).  We  can  conclude  A  —  AfSr  St  as  well. 

Hence,  letting  denote  submatrices  whose  value  is  irrelevant, 


Substituting  this  value  in  the  definition  of  Q'  given  in  Eq.  (11)  completes  the  proof: 


If  the  number  of  tokens  in  pi  is  always  at  least  m;  >  0,  we  can,  of  course,  define  S'  =  (m,-, . . . ,  nx-  —  1 } 
or,  equivalently,  change  the  definition  of  the  GSPN  so  that  the  range  of  tokens  in  p,  becomes  Sl  — 
{0, . . . , 7i j  —  77 ii  —  1}.  This  would  not  affect  the  proofs  in  this  paper,  but  could  improve  the  efficiency  of 
the  implementation. 
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Figure  2:  A  case  where  m^n,  n  6  <S,  but  m  <S. 

It  is  important  to  stress  that  and  U s^s^1  cannot  be  guaranteed  to  be  zero  (nor  can  R^^,, 

but  this  is  because  we  allow  timed  transitions  to  have  concession  in  vanishing  markings).  We  now 
formalize  this  observation,  already  implicit  in  previous  works  using  the  Kronecker  approach,  since  it  is 
fundamental  for  a  better  understanding  of  the  nature  of  the  matrices  R'  and  U'. 

Lemma  3.1  The  matrices  R'  and  TJ7  defined  in  Theorem  3.1  satisfy  the  following  “forward 
reachability  condition”: 

m  6  St  A  R'nn  >0  n  £  S  and  m  £  Sy  A  U'n  n  >  0  n  £  S  (15) 

However,  the  analogous  “backward  reachability  condition”  does  not  hold: 
n  £  «S  AR'n  ll  >0  7^  m  £  St  and  n  £  S  AUjn  u  >0  ^  m  £  Sy  (16) 


Proof:  It  is  straightforward  to  show  that  Eq.  (15)  holds  when  the  premises  of  Theorem  3.1 
are  satisfied.  We  simply  need  to  observe  that,  given  the  definition  of  R',  R ^  n  >  0  and 
m  £  St  imply  that  there  is  a  transition  tj  £  £(m)  and  that  its  firing  leads  to  n;  thus,  n  is 
reachable.  An  analogous  reasoning  holds  for  U'  and  m  £  Sy-  To  show  that  Eq.  (16)  holds,  it 
is  sufficient  to  give  an  example;  we  do  so  for  R'.  Consider  the  GSPN  in  Fig.  2,  having  positive 
finite  transition  rates.  The  firing  of  two  transitions  *3  could  lead  to  n  =  (0, 1, 1):  by  tf2, 
from  marking  (1,0,0),  and  by  ij,  from  marking  (1,0, 1).  In  our  notation,  rii  =  n2  =  n3  —  2, 


W1,1  = 


■  0 

0  ' 

1 

0 

W2’1  = 


'  0 

1  ■ 

0 

0 

w3’1  = 


■  1 

0  ‘ 

0 

1 

Thus,  the  contribution  of  U  to  R'101011  ls  K’^i.o’Woj  -W^J  =  w]  >  0,  hence  R'101  on  >  0. 
However,  given  the  initial  marking,  m°  =  (1,0,0),  the  marking  m  =  (1,0, 1)  is  not  reachable. 
□ 


Theorem  3.1  gives  a  characterization  of  the  infinitesimal  generator  Q  and  of  the  matrix  F  of  a  GSPN 
by  focusing  011  the  effect  of  each  transition  on  each  place.  An  alternative  statement  of  this  theorem,  is, 
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Figure  3:  A  GSPN  where  |<S7>|  <C  |<S'|,  and  one  where  |<S|  =  |«S' 


of  course, 


Q 

F 


+  RSt,Sv  (l|5u|  -  r'sl,Sv  ■  Usv,Si,)  U5y,5r  ~  ^ST,£ 
(l|Sv|  —  r’sv  ,SV  ‘ 


In  either  form,  however,  this  result  has  little  practical  value  in  itself,  since  both  expressions  contain 
an  inverse  which  cannot  be  expressed  using  Kronecker  operators  on  smaller  matrices.  One  case  where 
Theorem  3.1  has  a  direct  application  is  when  there  are  no  immediate  transitions.  Then,  St  =  S,  Sv  =  0, 
and  Eq.  (2)  simplifies  to  Q  =  R  -  A  =  RT>T  -  A  =  R^  5t  -  A'St  St. 

However,  a  solution  approach  based  on  this  idea  alone  has  limitations,  due  to  the  restrictions  that 
the  GSPN  must  satisfy.  Even  more  importantly,  though,  the  size  of  Q'  is  enormous,  potentially  leading 
to  inefficiencies.  Consider  for  example  the  GSPN  in  Figure  3(a),  having  positive  finite  transitions  rates. 
If  the  initial  marking  contains  a  total  of  n  tokens, 


|S|=(”  +  n  ')  «  |S’|  =  (»+1)‘. 

From  the  simple  case  when  n  =  1,  it  is  apparent  that  the  difference,  k  vs.  2k,  can  be  enormous.  For 
this  type  of  closed  networks,  Buchholz  [4]  suggested  a  solution  method  based  on  Kronecker  algebra 
that  does  not  create  unreachable  states,  applicable  when  the  interaction  between  submodels  is  of  the 
asynchronous  type  described  in  the  introduction. 

On  the  other  hand,  it  is  possible  for  S  to  equal  S'.  This  happens,  for  example,  in  a  live  free- 
choice  GSPN  with  capacities  whose  undirected  graph  obtained  by  ignoring  arc  directions  is  acyclic 
(this  is  a  generalization  of  [13,  Property  3],  which  refers,  however,  to  unbounded  marked  graphs). 
Another  example  is  that  of  open  acyclic  queueing  networks  with  communication  blocking  due  to  bounded 
buffers  [17]  which  could  be  named  “open  state  machines  with  capacities”  in  Petri  net  terminology.  The 
transitions  in  these  nets  have  at  most  one  input  and  one  output  place  and,  if  capacities  were  removed, 
every  place  would  become  unbounded.  See  the  GSPN  in  Figure  3(b)  for  a  simple  example  of  a  tandem 
network.  Indeed,  when  S'  =  St,  ’Fr(m)  is  simply  the  mixed-base  value  of  m,  hence  ^^(Jb)  does  not 
have  to  be  stored  explicitly.  Unfortunately,  such  a  situation  is  rare. 
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3.2  Merging  places  into  macroplaces 

The  type  of  marking  dependence  expressed  by  Eq.  (6)  and  (7)  is  quite  general,  but  for  example,  it  does 
not  let  us  specify  a  firing  rate  proportional  to  a  nonlinear  function  of  several  places  (e.g.,  minfm!,  m2}). 
We  now  show  how  this  limitation  can  be  overcome  in  practice  by  merging  places  (p,  and  p2,  in  our 
example). 

Consider  a  GSPN  A  =  (V,T  ,1,C~  ,C+ ,  m°,  w)  with  finite  reachability  set  S,  and  partition  V  into 
where  V,  =  {p,,,. .  .pijp  Then,  define  an  order-preserving  bijection  7  :  S  — »  S  C 

I 

7(mi,...,m|P|)  =  (mi,...,ih|^)  satisfying  7(111)  >-  7(11)  «=>  m  >-  n,  (17) 

where  m,  is  the  position,  in  lexicographic  order,  of  (m,-,, . . .  ()  in  the  set  obtained  by  projecting  S 

over  Vi. 

Lemma  3.2  Given  A,  V,  and  7  defined  as  above,  consider  the  “compacted”  GSPN 

A  =  (V,  t  =  T,  X  =  J,  C- ,  C+,  m°  =  7(m°),  w), 

where  the  input  and  output  arc  cardinalities  are  defined  to  ensure  that,  in  corresponding 
markings  m  and  7(m)  =  m,  tj  €  T  has  concession  in  A  iff  it  has  concession  in  A  and  that, 

in  this  case,  m-^n  =  7(11)  in  A  iff  m-^n  in  A  ,  while  the  weights  for  tj  are  defined  to  have 
the  same  value  in  corresponding  markings: 

•  If  tj  €  £(m)  and  its  firing  does  not  change  the  marking  of  any  place  in  Vi,  that  is,  if 
Vpi  €  Vi  :  C,"  (m)  <  m,  A  C,"  (m)  =  C,+ (m),  define  C“;(m)  =  C+  (m)  =  0. 

•  If  tj  e  £(m)  and  its  firing  changes  the  marking  of  some  place(s)  in  Vi,  that  is,  if 
Vpi  €  Vi  :  Cfj(m)  <  m/  A  3p/  €  pi,  Cz~  (m)  /  Cj^(m),  define  C“y(rh)  =  m;  and 
Ct.(m)  =  hi. 

•  Otherwise,  tj  is  disabled  in  m,  that  is,  3p/  €  Vi,  C(“  ( m)  >  m(;  then,  define  C^m)  = 
m,  +  1,  while  the  value  of  C^(m)  is  irrelevant. 

•  Define  Wj(m)  =  Wj(m). 

Then,  the  transition  rate  matrices  R  and  R,  defined  by  A  and  A,  respectively,  are  identical. 

Proof:  Omitted  for  brevity  (it  is  sufficient  to  show  that  the  stochastic  processes  described 
by  A  and  A  are  identical).  □ 

Lemma  3.2  allows  us  to  compact  an  arbitrary  set  of  places  into  a  single  place,  which,  together 
with  the  transitions  connected  to  it,  corresponds  to  a  sub-GSPN  of  [11,  12].  This  operation  must  be 
performed  when  the  marking  dependencies  in  the  GSPN  are  not  of  the  type  allowed  by  Theorem  3.1. 
It  might  be  performed,  even  when  the  theorem  is  applicable,  to  reduce  the  number  of  matrices  involved 
in  the  description  of  R  at  the  cost  of  increasing  their  size. 

From  now  on,  macroplaces  are  indicated  as  dashed  boxes  surrounding  sets  of  places;  the  compacted 
GSPNs  are  not  shown  explicitly,  since  they  would  not  add  to  the  comprehension  of  the  model. 
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4  Timed  synchronizing  transitions 

The  contributions  in  [4,  5,  11,  12]  have  assumed  that  the  GSPN  is  decomposed  in  such  a  way  that  each 
iteration  of  the  solution  method  performs  Kronecker  products  of  few  large  (but  manageable)  matrices, 
while  Theorem  3.1  uses  many  (|T|  -  \P\ )  small  (n;  x  n:)  matrices. 

Lemma  3.2  addresses  the  size  issue:  we  can  merge  places,  thus  obtaining  \T\  •  \V\  larger  matrices.  In 
this  section,  we  show  how  the  number  of  matrices  involved  can  be  further  reduced  by  merging  transitions, 
or  rather,  the  corresponding  matrices.  The  results  are  similar  to  those  derived  by  previous  authors,  who 
assumed  all  synchronizing  transitions  are  timed,  but  we  present  them  here  for  three  reasons.  First, 
we  exhibit  substantially  different  proofs  for  these  results;  [4,  5,  11,  12]  consider  a  set  of  sub-GSPNs 
and  combine  them  using  synchronizing  transitions,  thus  Kronecker  operators  are  introduced  only  at 
the  last  step.  Instead,  we  start  from  the  Kronecker  expression  of  Theorem  3.1  for  the  entire  GSPN 
and  derive  our  results  by  exploiting  the  properties  of  Kronecker  operators.  Second,  our  results  include 
the  management  of  immediate  transitions  and  vanishing  markings,  while  previous  works  have  simply 
assumed  that  these  are  eliminated  locally  using  the  traditional  approach.  Finally  and  most  importantly, 
our  result  apply  to  a  larger  class  of  GSPNs,  since  a  more  general  marking-dependent  behavior  is  allowed 
by  Theorem  3.1. 

4.1  Partitioning  the  set  of  transitions 

Without  loss  of  generality,  we  assume  from  now  on  that  each  transition  in  T  has  at  least  one  input  or 
one  output  place:  Vtj  €  T,3pt-  €  V,  ^  0  V  ^  0.  Then,  let  T*  C  T  be  the  set  of  “local” 
transitions  which  affect,  or  are  affected  by,  only  a  single  place  p,: 

Vpi  €  V  :  T*  =  {tj  €  T  I  Vp,  e  V,pi  f  p„  C,“.  =  C+  =  0  A  u,()i  =  l}  ,  (18) 

and  T*  =  T\  UPlgp  T*  be  the  set  of  synchronizing  transitions  which  instead  affect  or  are  affected  by  at 
least  two  places.  Clearly,  these  sets  constitute  a  partition  of  T.  Also,  let  X*  =  T*  n  X ,  J*  =  T*  HI, 
X*  =  T*  fl  X,  and  T  =  Ti  D X. 

Lemma  4.1  Consider  a  GSPN  satisfying  the  requirements  of  Theorem  3.1.  Then, 


R'=  £ 

(g)  vr* 

+  R’  , 

u'=  £  »; 

•  (g  w*> 

*+  ©u\ 

(19) 

Pier 

Pi£V 

Pi€T> 

Pi&'p 

A'=  T,  wj 

•  (g)  A1'1 

+  ©  A'  , 

V 

w 

II 

•  <g>  riJ 

+  ©  r\ 

(20) 

Pi€V 

Pi€P 

ij€l* 

ViZV 

where  R! ,  U‘,  A\  and  are  square  matrices  of  size  n,  x  n,  defined  as 

r>  =  £  w; .  wiJ  ,  u*  =  ,  x  =  Y,  wr  A'J  .  T=  £ 

i,ei<  tje.v  tjep 

Hence,  as  special  cases,  R‘  =  A'  =  0  if  X ’  =  0  and  U‘  =  J"  =  0  if  2‘  =  0. 
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Proof:  We  only  prove  the  result  for  R',  the  proof  for  the  other  matrices  is  analogous.  Given 
the  condition  specified  by  Eq.  (18)5  we  know  that  W/j  =  ini  if  p{  p{  and  tj  €  T\  Then, 
the  proof  is  a  simple  matter  of  matrix  manipulation  within  the  Kronecker  expressions: 


R' 


E  wr 

0  WiJ 

t3ZX 

Pi€T> 

E 

wi 

•  0  W-J 

f  +  E  E  ws 

■  0  WlJ 

‘>e-v 

p.ev 

Pi€'PtjSXi 

piev 

E 

w* 

•  0 

f+  E  e*; 

■  In,  ®  • 

•  •  ®  Inj_,  <8>  W’  J  <g)  Ini+]  (g>  •  -  •  ®  Injpj 

tj€X* 

P.eP 

Piev  tjex* 

/  \ 

E 

w] 

•  0  W*-J 

®  [  E  «$-Ww)  <S»I„i+1 

P.e? 

P.6J' 

V^e.v*  / 

E 

w* 

•  0  wij 

;+  E1..®- 

®  I»;  _  ] 

®  R1  ®  In,+1  ®  ®  Iti|P| 

P.e? 

PiPV 

E 

WJ 

•  0  W<lJ 

1  +  ©  R' 

tyeA'* 

p.ev 

p.ev 

This  partition  reduces  the  number  of  Kronecker  product  terms  from  \X\  to  \X*\  in  R'  and  from  \1\  to 
|J*|  in  U',  respectively,  and  adds  one  Kronecker  sum  to  both.  Transitions  satisfying  Eq.  (18)  arise  after 
applying  Lemma  3.2,  that  is,  after  “decomposing”  a  large  GSPN  into  several  smaller  sub-GSPNs.  Each 
sub-GSPN  corresponds,  in  our  terminology,  to  a  (macro)place,  plus  the  set  of  transitions  local  to  it. 
This  transformation  does  not  have  to  be  explicitly  performed  in  practice,  only  its  result,  the  matrices 
corresponding  to  the  set  of  macroplaces,  need  to  be  computed.  A  good  partition  of  the  places  results 
in  a  compacted  GSPN  where  most  transitions  are  local  and  the  number  of  tokens  in  each  compacted 
place  (number  of  markings  in  the  sub-GSPN,  in  the  terminology  of  [11,  12]),  is  manageable.  Methods 
to  determine  a  good  partition  are  beyond  the  scope  of  this  paper  and  are  left  for  future  research. 


4*2  An  efficient  Kronecker  expression  for  the  CTMC 

Given  any  GSPN,  we  can  always  apply  Lemma  3.2,  resulting  in  a  compacted  GSPN  satisfying  the 
requirements  of  Theorem  3.1,  then  apply  Lemma  4.1.  If  the  partition  is  such  that  all  immediate 
transitions  are  local,  we  can  then  restate  the  main  results  of  [11,  12]  in  a  more  general  setting. 

Theorem  4.1  Consider  a  GSPN  satisfying  the  same  requirements  as  for  Theorem  3.1,  and 
such  that  all  immediate  transitions  are  local: 

r-l  =>  U^QU1,  T'  =  X\  (21) 

pi€V 


Then 


Q"  =  X>;  •  (g>  (W-  •  X‘)  +  0  (R*  •  X’)  -  A', 

ij€X‘  pier  Pi€V 
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where 


X*  =  (lni  -  (T*  +  ry'lf)  1  r  =  rowsum(U‘)  T*  =  In.  -  S(r), 

satisfies  Q ^  5  =  Q,  with  the  meaning  defined  in  Eq.  (2). 

Proof:  First,  we  observe  that  condition  Eq.  (21)  implies  J1'  =  ®p.€7>  -T4  and  T'  =  0Piej>  T4. 
In  other  words,  a  “global”  marking  m  is  tangible  iff  all  its  “local”  components  are  tangible. 
We  can  then  manipulate  Q"  as  follows: 

Q"  =  £  •  0  (w4li  •  X‘)  +  0  (R’  •  X4)  -  a! 

tj&X'  Piev  Pier 

=  E  ®w“.  ®  X'  +  ®  R'.  ®  r  +  $  (tf  r)  -  $R‘. 

tj€X'  Pi€T  pit?  P,er  Pier  Pier  Pier  Pier 

' - v - ' 

D' 

£  w*  •  (g)  Wij'  +  0  R*  j  •  (g)  X‘  +  D'  -  A! 
hex*  Pier  p&r  )  Pier 

=  R'  •  (g)  X4  +  D'  -  A! 

Pier 

Partition  S'  into  <Sf  and  <S(/,  corresponding  to  local  markings  enabling  only  timed  local 
transitions,  or  some  immediate  local  transition,  respectively,  and  rearrange  the  rows  and 
columns  of  X'  accordingly: 

x‘  =  (in,  -  (t*  +  r4)-1u4')_1  = 

where  the  subscripts  “T,  T”,  “V,T”,  and  “V,  V”  have  the  usual  meaning,  but  applied  to  the 
local  matrix  for  place  i.  =  (i^,^  —  (r'sv,sv)~l  describes  the  expected 

number  of  visits  to  each  local  vanishing  marking  before  reaching  a  local  tangible  marking, 
starting  from  each  each  local  vanishing  marking,  while  Ply>T  =  •  (riSv  Sv)~'  •  s 

describes  the  probability  of  reaching  each  local  tangible  marking,  starting  from  each  local 
vanishing  marking. 

We  continue  assuming  that  \V\  =  2,  the  general  proof  follows  exactly  the  same  idea.  Local 
matrices  are  partitioned  according  to  whether  the  corresponding  local  markings  are  tangi¬ 
ble  or  vanishing  (regardless  of  whether  the  global  markings  are  reachable  or  not).  Global 
matrices  are  partitioned  according  to  the  following  order:  tangible  states,  vanishing  states 
enabling  only  immediate  transitions  in  T2,  vanishing  states  enabling  only  immediate  transi¬ 
tions  in  T1 ,  and  vanishing  states  enabling  immediate  transitions  in  both  T1  and  T2.  First, 
we  show  that  the  tangible  rows  of  D'  are  zero: 

D'  =  (R’X1  ©  R2X2)  —  (R1  ©  R2)(X’  ®  X2) 

=  R'X1  0  In2  +  In,  <g>  R2X2  —  R’X1  ®  X2  —  X1  ®  R2X2 


Tl 

aT,T 

0 

p  i 

.  rV.T 

K.v . 
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Figure  4:  A  GSPN  lo  illustrate  the  difference  between  Q'  and  Q". 
=  R’X1  ®  (I)l2  -  X2)  +  (!»,  -X,)®R2X2 


since  X1  and  Ini  have  the  same  tangible  rows.  Then,  we  show  that  the  tangible  columns  of 
®Vi€P  X*  and  (I  -  (T'  +  r'y'U'y1  coincide: 


The  blocks  in  the  first  (tangible)  column  of  this  last  matrix  contain  the  correct  values, 
since  the  top  left  block  is  simply  the  identity  and  the  other  blocks  correctly  describe  the 
probabilities  of  reaching  tangible  markings  from  vanishing  markings,  which  are  the  values 
in  the  corresponding  blocks  of  (I  —  (T'  +  -T')-1!!')-1 .  From  Eq.  (22)  and  (23)  we  can  then 
conclude  that  Q$T  $T  =  Q,  as  in  the  proof  of  Theorem  3.1.  □. 

In  practice,  Theorem  4.1  is  used  to  generate  only  the  relevant  portion  of  Q"  in  the  numerical  solution 
method.  In  other  words,  we  eliminate  the  vanishing  markings  “on  the  fly”  ( as  in  [4,  5,  10,  11,  12]) 
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We  observe  that  both  Q'  and  Q"  describe  Q,  but  the  two  normally  differ.  For  example,  consider  the 
GSPN  in  Fig.  4,  having  a  single  tangible  state  and  three  vanishing  states  (the  example  is  trivial,  but  it 
is  sufficient  to  illustrate  the  point).  Assuming  that  the  rates  of  timed  transitions  ^3,  and  is  are  a,  6, 
and  c,  respectively,  and  that  the  weights  of  transitions  t2  and  1 4  are  2  and  3: 


R1  = 


■  0  «  ' 
0  0 

R2  = 

0  b 

0  0 

w1-5  = 

or 
0  0 

w2'5  = 

or 
0  0 

u1  = 

0  0' 
2  0 

u2  = 

0  0 

O  CO 

■  ■ 

From  which  we  obtain 


W1  = 


'00' 
0  2 

w2  = 

00' 

0  3 

T1  = 

10' 
0  0 

(j»2  _ 

'  1  o' 
0  0 

X1  = 

'10' 

11 

X2  = 

10' 

11 

and 


R' 


"  0 

b 

a 

c 

■  1 

0 

0 

0  ■ 

'  0 

0 

0 

0  ' 

'  0 

0 

0 

0  * 

0 

0 

0 

a 

0 

0 

0 

0 

U'  = 

3 

0 

0 

0 

y  1  / 

0 

3 

0 

0 

T'  = 

0 

0 

r  - 

0 

0 

0 

0 

0 

b 

0 

0 

0 

0 

2 

0 

2 

0 

.  0 

0 

0 

0  . 

.  0 

0 

0 

0 . 

.  0 

2 

3 

0  . 

.  0 

0 

0 

5 . 

The  resulting  Q'  and  Q"  matrices  are  then  different: 


Q'  = 


a  +  b+  c 

b  +  \c 

a  +  |c 

c  " 

a  +  b  +  c 

b  +  c 

a  +  c  c 

a 

\a 

fa 

a 

0 

a 

0  a 

b 

5 

¥ 

5 

¥ 

b 

-A! 

Q"  = 

0 

0 

b  b 

0 

0 

0 

0  . 

0 

0 

— 1 

0 

0 

—  A'. 


Indeed,  the  difference  between  Q'  and  Q"  is  already  apparent  from  Eq.  (23).  The  diagonal  blocks 
Ij  T  ®  y  and  N yV  ®  I^  T  correctly  describe  the  expected  number  of  times  the  corresponding  global 
vanishing  markings  (enabling  only  local  immediate  transitions  in  T2  or  T1,  respectively)  are  entered, 
given  that  a  timed  transition  firing  leads  to  the  corresponding  diagonal  block.  However,  the  last  three 
blocks  on  the  bottom  row  do  not  reflect  the  same  quantities  when  a  timed  transition  firing  leads  to 
vanishing  markings  enabling  immediate  transitions  in  both  T1  and  T2.  In  particular,  P{,T  ® 
describes  the  correct  quantity  only  if  we  could  assume  that  all  enabled  immediate  transitions  in  T1  keep 
firing  before  any  of  those  in  T2  do,  which  is  not  necessarily  the  case,  while  Nj,  v  ®  PyT  assumes  the 
opposite.  Finally,  N]Vy  <g>  N \v  does  not  reflect  the  number  of  times  global  markings  are  entered  at  all. 
This  leads  us  to  the  following  observation. 


Corollary  4.1  If  the  GPSN  of  Theorem  4.1  is  such  that  the  firing  of  any  timed  (synchro¬ 
nizing)  transition  tj  6  A'#  enables  (local)  immediate  transitions  in  at  most  one  set  Tl,  for 
some  pi  £  P,  then  Q sTtsv  =  F,  as  defined  in  Eq.  (5). 

Proof:  The  condition  for  this  corollary  implies 

C  U  S±  x  •  •  •  x  4-1  x  Si  x  4+1  x  •  •  •  x  s?1. 

Pi€ V 
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As  discussed  in  Theorem  4.1,  the  blocks  on  the  columns  corresponding  to  this  type  of  vanish¬ 
ing  markings  are  computed  correctly  in  see  Eq.  (23).  The  bottom  rows  of  Eq.  23, 

corresponding  to  any  (unreachable)  vanishing  markings  m  enabling  immediate  transitions 
in  more  than  one  set  T*,  are  irrelevant,  since  Lemma  3.1  guarantees  that  R^,  m  =  0  in  this 
case.  D 


5  Immediate  synchronizing  transitions 

If  some  of  the  synchronizing  transitions  are  immediate,  Theorem  3.1  still  applies,  but  Theorem  4.1, 
which  allows  the  efficient  computation  of  the  solution  in  practice,  does  not.  In  this  section,  we  show 
how  “preservation  of  the  vanishing  markings”  [8]  can  be  used  to  remove  this  limitation,  also  present  in 
[11,12]. 


5.1  Embedding  a  DTMC 

First,  we  summarize  the  main  ideas  in  [8],  which  examines  an  alternate  method  to  compute  7r: 

•  Define  the  transition  probability  matrix  P  of  the  embedded  DTMC,  expressing  the  probability  of 
going,  in  one  firing,  from  any  marking  m  €  S  to  any  other  marking  n  e  S,  regardless  of  whether 
they  are  tangible  or  vanishing: 

p  _  A*  Rr.r  I  1  (24) 

\JVyT  VvyV 


Compute  the  steady-state  probability  vector  7  €  JR)S^  of  the  embedded  DTMC: 

7  •  P  =  7  subject  to  the  normalization  7  •  l|s|xi  =  1- 


Obtain  both  7r  and  <f>  from  7,  using  the  holding  times  in  the  tangible  markings  as  weights: 


Vm  €  St,  7rm  = 


7m  '  hm 
£neST  T11  •  hn 


and  Vm  G  5,  =  — - 111 — — .  (26) 

2^ii65t  7n  '  “n 


The  correctness  of  the  method  can  be  verified  by  observing  that,  from 

r  1  1  r  A_1Rt,v  1  r  I  1 

[It  I  7v]  •  — 77 - 77 -  =  [7 T  I  7v]> 

\JV,T  Uv,V 


we  can  obtain  7 v  =  7Tyl  (I  —  Uyy)  1 ,  and,  substituting  it  in  the  above  equation, 

7r/l_1  •  (R7yr  +  R7’,v  (I  —  U^r1  U v/y)  =  0. 

" - ^  "  .  ' - - - - 

TT-^-yi  l  i|i'|Xij  Q 
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We  can  then  divide  both  sides  of  the  equation  by  the  constant  7T  •  A~l  •  l|s|xi,  resulting  in  tt  •  Q  =  0. 

In  [8],  it  was  found  that  the  solution  time  is  often  greater  than  with  the  “elimination”  approach 
based  on  Eq.  (3).  This  is  due  to  the  number  of  nonzero  entries  in  P,  normally  larger  than  in  Q,  and, 
frequently,  to  a  slower  numerical  convergence.  However,  pathological  cases  where  P  has  substantially 
fewer  entries  than  Q  arise  when  N  tangible  markings  can  reach  a  small  set  of  vanishing  markings,  which 
can,  in  turn  reach  M  tangible  markings.  This  uN-to-M  switch”  behavior  corresponds  to  0(N  +  M) 
arcs  in  P  and  0(N  ■  M)  in  Q,  hence  it  affects  the  elimination  approach  negatively  both  in  terms  of 
storage  and  execution  time,  although  the  number  of  iterations  in  the  numerical  solution  might  still  be 
smaller  with  elimination.  The  use  of  preservation  results  in  the  following  analog  of  Theorem  3.1. 

Theorem  5.1  Under  the  same  conditions  of  Theorem  3.1,  the  matrix 

P'  =  (T'  •  A'  +  r')-1  •  (T'  •  R'  +  U')  (27) 


satisfies  P^5  =  P,  as  defined  in  Eq.  (24). 

Proof:  The  pre-multiplication  of  A'  and  R'  by  T'  eliminates  the  effect  of  timed  transitions 
having  concession  in  vanishing  markings.  However,  if  we  focus  on  the  reachable  states,  the 
statement  of  the  theorem  is  equivalent  to  saying  that 


P  = 


A'< 


-l 


S'p^S'p 


■  R'< 


=T 


‘S'p  ,5 


r'< 


Sy,Sv 


•U  'SvtS 


(28) 


This  equality  then  follows  from  the  definition  of  P  and  the  meaning  of  R'  and  U'  already 
established  in  Theorem  3.1.  Eq.  28  is,  of  course,  the  expression  used  in  practice  for  a 
numerical  solution.  D 


The  efficiency  of  a  solution  based  on  Eq.  (28)  is  improved  by  exploiting  the  existence  of  local 
transitions  (Lemma  4.1): 


P  = 


(  £  ®  A<,i  +  ©  a1)  •  (  £  “J-  <8> W<J  +  ©  R1 

psp  ner  /  V,ca-  ner  ner  ) 

_  _i  7  ; 


ST:S 


E  •  <g>  +  ©  r 

[  P.6P  p.€P  ) 


•  E  -0^+0  ir 

Vb6l‘  Piev  vi€V  ) 


Sv,S 


(29) 


since  this  expression  for  P  reduces  the  number  of  Kronecker  products  to  be  performed  at  each  iteration. 


5.2  Using  partial  elimination  to  improve  solution  efficiency 

A  disadvantage  of  the  approach  just  described  is  that  the  size  of  the  probability  vector  7  for  the  DTMC 
is  now  |<S|,  considerably  larger  than  |5r|  in  many  practical  models.  Analogously,  the  size  of  the  matrices 
for  place  i  is  given  by  the  projection  of  S  onto  its  i-th  component,  <S‘  =  {/  :  3m  €  S,  m,  =  /},  regardless 
of  whether  markings  satisfying  m,  =  l  are  vanishing  or  tangible. 
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It  was  already  observed  in  [8]  that  it  is  possible  to  eliminate  a  subset  of  the  vanishing  markings, 
preserving  only  those  involved  in  large  switches,  in  the  hope  of  achieving  the  best  memory-execution 
tradeoff.  We  can  exploit  the  same  idea  in  our  approach,  but  for  a  different  purpose.  Partition  the  local 
state  space  for  place  i,  S' ,  into 

•  =  {l :  Vm  €  St,  m ,•  =  /},  the  set  of  local  tangible  markings. 

•  S's  -  {l  :  3m  e  Sv,mi  =  l  A  T  G  £(m)  ^  0),  the  set  of  possibly  synchronized  local  vanishing 
markings. 

•  S[  =  {l  :  Vm  €  «Sv,m,  =  /  =>  I*  D  £(m)  =  0},  the  set  of  non-synchronized  local  vanishing 
markings. 

Any  immediate  synchronizing  transition  can  be  enabled  only  in  markings  having  components  in  U <S^ . 
that  is,  in  Ss  =  <Sv  H  ((<Sj  U  5$)  x  •  •  •  x  (S|^  U  £5**)) . 

We  now  define  a  “partially  eliminated”  (or  “partially  preserved”)  DTMC  with  transition  probability 
matrix  P  and  state  space  Sk  =  StUSs,  (I<  stands  for  ‘keep”)  which  can  be  used  to  compute  7r  and  4> 
more  efficiently  than  from  P.  Partition  the  matrices  Wij‘,  R‘,  U‘,  A''3,  T'  3 ,  A',  and  I”,  according  to 
the  sets  S'K  =  S^  U  S's  and  S[.  For  example, 


Then,  assuming  that  the  rows  U‘LiA-|U‘LjLJ  are  already  normalized  (this  can  be  easily  enforced  since 
each  U*  is  built  before  starting  the  overall  solution),  define  the  matrices 

W*  =  w%K  +  W&  •  (I  -  u^) •  u^, 

R'  =  R‘k  k  +  R 'K,L  ■  (l  ~  •  U^,and 

U‘  =  u^  +  ujf(L.(i-ui>t)"1.ui,K. 

Given  this  definition,  the  blocks  for  the  rows  and  columns  of  Sk  in  A''3 ,  Jn!J,  A\  and  P*  still  contain 
the  correct  row  sums  for  the  corresponding  matrices.  Then,  we  can  state  our  final  theorem,  which 
allows  the  efficient  solution  of  a  structured  GSPN  with  immediate  synchronizing  transitions. 

Theorem  5.2  Under  the  same  conditions  of  Theorem  3.1,  define  the  transition  probability 
matrix 


P 

(30) 


E  w*j  (8)  AKtK  +  ®  a'k  k 

tjEX-  p.€P _ Pi  6  V  )_ 

E  •  <8>  rti<  +  ®  rKJ< 

Pi€V  Pi€? 


.*>€*•  vie?  pier  /  St,s,, 


E  <g>  wiJ+  0  u1' 

(,ei*  p,€V  Pi&v 
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and  solve  for  7  •  P  =  7  subject  to  the  normalization  7  •  l|sA.|xl  =  1.  Then,  the  steady-state 
probability  and  the  rate  of  entering  the  markings  in  Sk  are 


Vm  €  St,  7rm  = 


7m  '  h 


m 


UnCSy  7n  ^11 


Vm  e  SK,  4>m  = 


7m 


EneST  7n  ■  hn 


Proof:  We  only  need  to  show  that  P  correctly  describes  the  DTMC  obtained  when  embed¬ 
ding  the  GSPN  at  the  times  when  markings  in  StUSs,  but  not  those  in  S 1  =  S\(StUSs), 
are  entered.  In  other  words,  7  should  differ  from  -fK  only  by  a  multiplicative  constant, 
where  [7^  |  7^]  =  7  is  the  solution  to  Eq.  25,  which,  in  block  form,  is 


[7aI7l] 


P  A, A 

P  A,L 

P  L,K 

P  L,L 

[TaItJ- 


We  can  then  obtain  7^  =  7^  •  Patl  •  (I  —  Pl,l)  1  and,  by  substitution, 


7  a  •  (Pa'.l  •  (I  -  P l,l)  1  •  Pl.a  +  Pa, a)  —  7a- 

Then,  it  is  sufficient  to  show  that  (Pa.l  ■  (I  -  P l,a)-1  •  P l,k  +  Pa, a)  and  P  coincide.  For 
any  m,n  €  Sk,  (Pa,a)„m1  represents  the  probability  of  going  from  marking  m  6  Sk  to 
marking  n  €  Sk  in  a  single  firing,  while  (Pa.l  •  (I  -  Pa, a)-1  •  PA,A')m>n  represents  the 
probability  of  going  from  m  to  any  marking  m1  €  Si,  visiting  any  number  of  markings  in 
Si,  and  finally  leaving  Si  from  some  marking  m2  (possibly  the  same  as  m1)  to  reach  n  in 
one  firing. 

Assuming  that  m  £  Sr  and  n  differs  from  m  in  at  most  the  position  for  place  pi,  the 
corresponding  entry  Pm,n  is 


P  w,n  — 


=  A 


=  A 


(  <g>  Ati<  +  0  ^a,a)  •  (  £  ■  0  +  0  R’) 

\tj6.i'*  p.ev  Pier  j  m  m  \t,zx'  Piev  K&>  /m  n 

1,111  ‘  f  E  •  II  . 

\h  ex*  Pier  ) 

:,».•(  E  »;•  n  (wi*+  e  ^  e 

\be-v  pitr  \  m;€5[  m?e5[  ■ 

+  +  E  »!»„»,'  E 

V  ml*S,L  rnjesi  **  '  J  j 


(if  n  and  m  differ  in  more  than  one  position,  the  cause  must  be  the  firing  of  a  synchronizing 
transition,  so  the  “local  term”  for  place  p;  in  the  last  expression  is  absent).  In  any  case, 
Ahlm  is  just  a  normalization  factor  (if  m  £  Ss,  r~^m  would  be  used  instead),  so  the 
expression  indicates  the  required  probability.  The  key  issue  is  that  the  order  of  firing  of 
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diate  transitions.  For  example,  consider  the  GSPN  in  Fig.  5(a).  We  have  £?>  =  {(0),(2)},  S$  =  {(1)}, 
Si  =  0,  S'l  =  {(0)},  c>|  =  {(1)},  and  Si  —  {(2)},  as  shown  in  Fig.  5(b).  After  eliminating  the  local 
markings  Si  from  the  second  sub-GSPN,  we  obtain  the  “local  reachability  graphs”  in  Fig.  5(c).  Fig. 
5(d),  5(e),  and  5(f)  show  the  graphs  describing  P',  P,  and  P,  respectively.  The  dotted  arcs  in  Fig.  5(d) 
correspond  to  timed  transitions  with  concession  in  vanishing  markings,  which  lead  to  markings  in  S'. 
These  are  not  present  in  P'.  Hence,  (global)  marking  (1,2)  is  unreachable  and  is  absent  in  P  and  P. 
Furthermore,  markings  (0,2)  and  (2,2)  are  absent  from  P,  because  their  second  component,  2,  enables 
only  local  immediate  transitions  in  the  second  sub-GSPN,  (2)  £  «S§.  However,  marking  (1,0)  is  still 
present  in  P,  even  if  it  enables  only  t3,  a  local  immediate  transition.  This  is  because  its  first  compo¬ 
nent,  1,  corresponds  to  having  a  token  in  p3,  which  is  a  condition  for  the  enabling  of  the  synchronizing 
immediate  transition  t7.  In  other  words,  we  cannot  eliminate  the  local  marking  (1)  from  the  local  state 
space  for  the  first  sub-GSPN,  because  this  would  eliminate  both  global  markings  (1,0)  and  (1, 1),  and 
eliminating  (1,1)  would  make  it  impossible  to  capture  the  efTect  of  synchronizing  transition  <7  in  the 
Kronecker  products  of  Eq.  (30). 

6  Numerical  solution 

Applying  Theorem  4.1,  we  use  the  generator  matrix  Q  =  QsTist  to  compute  the  stationary  distribution 
7r  6  IR}St\  of  the  CTMC  underlying  the  GSPN  according  to  Eq.  (3).  We  use  an  approach  based 
on  Kronecker  algebra  to  avoid  storing  Q  explicitly,  so  only  iterative  methods  which  do  not  require  the 
modification  of  Q  itself  can  be  used  effectively.  Adopting  the  Jacobi  method  with  overrelaxation  (JOR), 
we  transform  Q  into  the  iteration  matrix  M  =  (1  —  cu)  •  I  +  u>  •  R  •  diag(h)  and  solve  the  eigenvector 
problem: 

7T  -  M  =  7T  subject  to  7T  •  1|5T|X1  =  1. 

Successive  approximations  of  7r  are  obtained  iteratively  as 

Tr^+il  «_  tj-M  .  M,  (31) 

starting  from  an  initial  probability  vector  7rl°l  satisfying  7rl°l  >  0  and  7r^l|5T|Xi  =  1.  If  the  CTMC  is 
ergoclic  and  if  the  iterations  converge,  JOR  is  guaranteed  to  result  in  the  correct  solution,  regardless 
of  the  value  of  tt^,  that  is,  i r  =  lirnm_*oo^m^  if  this  limitexists.  We  do  not  discuss  the  choice  of  the 
relaxation  parameter  u;,  0  <  w  <  2,  which  affects  the  convergence  rate  [16];  for  a  detailed  analysis  of 
numerical  techniques  for  the  solution  of  Markov  chains,  see  [15]. 

We  then  need  to  multiply  a  full  vector,  by  a  matrix,  R,  which  is  described  as  a  (submatrix  of 
a)  Kronecker  expression  of  smaller  matrices, 

r=  ( E  igiw^+eir] 

Pi6P  /  St  St 

where  R!  =  Rl  •  X*  and  W’,J  =  W1,J  •  X‘.  A  similar  discussion  applies  when  using  Corollary  4.1  (to 
obtain  <pnl  for  m  €  Sy  after  7r  has  been  obtained,  we  compute  7r  •  QsTtsv )  or  Theorem  5.2  (to  obtain 
7^,n+1l,  we  perform  the  product  7^  •  P). 
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•  During  the  reachability  set  construction,  St  is  normally  stored  in  lexicographic  order  in  a  balanced 
search  tree.  However,  we  can  assume  from  now  on  that  is  accessed  exclusively  using  'Pj1. 

•  After  the  reachability  graph  construction,  is  stored  as  an  array  of  size  \St\  of  pointers  to  the 

markings,  points  the  k- th  marking,  that  is,  the  marking  m  satisfying  'l'(m)  =  k.  Several 

methods  can  be  used  to  store  a  marking.  In  principle,  a  single  integer  describing  the  position  of 
m  in  S'  according  to  lexicographic  order  is  sufficient,  but,  in  practice,  this  integer  might  require 
more  than  32  bits. 

•  h  €  carries  the  same  information  as  A.  Instead  of  storing  this  vector  explicitly,  we  could 

recompute  A  at  each  iteration,  to  further  reduce  the  memory  requirements.  This  is  a  memory- 
execution  trade-off. 


It  is  important  to  note  that  an  alternative  method  suggested  by  Kemper  [12]  for  the  state  space  ex¬ 
ploration  has  the  advantage  of  allowing  the  determination  of  whether  a  marking  is  reachable  in  0(1) 
instead  of  0(log  \St\)  operations.  However,  it  requires  a  bit  vector  of  size  |5'|,  which  might  be  a  problem 
when  \S'\  >  |«ST|. 

If  we  define 

i- i  PI 


Vpi  6  V ,  low(i )  =  II and  up(i )  =  JJ  n/, 


<=i 


l=i+ 1 


the  contribution  of  the  Kronecker  sum  to  the  new  iteration  vector  7r^m+1J  in  Eq.  (31)  can  be  rewritten 
as 


tth  •  0  R*  =  E  *h  •  (W)  ®  Qi  ®  i*p(o) 

Pi£V  p;€P 


E 

P.eP 


while  the  contribution  of  the  Kronecker  products  is 


E 

be  a- 


IV- 


7Tl 


0  w1 
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TV 


[m] 


PI  , 

n(i 

i-l 
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®  w 


®  ^(o)  =  E 

qe.v« 


W; 


7r[,n+11. 


We  implemented  the  basic  operation  a  *  (I/0™(i)  ®  A  <g>  Iup(i))  in  a  function  mulifa  A,i),  where  a  is 
a  vector  of  size  |5r|  and  A  is  a  rii  x  n,  matrix.  Each  term  is  computed  with  a  single  call  to 

mult ,  while  each  term  7r[m+^  requires  \V\  successive  calls  to  mult .  During  this  iteration,  some  states 
(possibly  unreachable)  are  passed  and  temporarily  stored.  However,  if  all  the  states  of  S'  are  reachable 
and  tangible,  that  is,  if  [^S'l  =  |<St|,  there  is  no  need  to  explicitly  identify  the  reachable  states.  We  stress 
that  this  definition  allows  any  number  of  Kronecker  products  to  be  performed,  so  it  does  not  restrict 
the  possible  number  of  sub-GSPNs  involved  in  a  synchronization. 

For  the  numerical  computation  we  consider  a  fork /join  kanban  network  of  four  sub-GSPNs  as  given 
in  Figure  6(b).  Each  sub-GSPN  i  =  1,...,4,  of  the  network  is  modelled  by  the  GSPN  shown  in  Fig. 
6(a)  on  the  left.  Multiplying  R*  and  W*’-7  by  X*  corresponds,  de  facto,  to  the  automatic  elimination 
of  the  immediate  transitions  tredoi  and  tQk resulting  in  the  GSPN  shown  in  Fig.  6(a)  on  the  right, 
where  the  rates  of  tmiredoi  and  im|.0*t.  are  the  the  rate  of  tm  multiplied  by  the  firing  probabilities  of  tredo , 
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Figure  6:  A  sub-GSPN  and  a  fork/join  arrangement  of  four  sub-GSPNs. 

and  toki,  respectively.  The  parameters  specifying  the  stochastic  behavior  of  the  sub-GSPNs  are  the 
rates  w,nj,  w„w,  and  W(,ac^  and  the  probabilities  w0/i.i  and  w^.  We  assume  that  these  quantities  are 
constant,  but  they  could  depend  on  the  local  marking  (or  even  on  the  global  marking,  in  the  case  of 
synchronizing  transitions)  without  affecting  the  feasibility  or  the  complexity  of  the  solution  algorithms. 
Pallets  enter  sub-GSPN  i  of  the  kanban  network  through  transition  tinn  which  requires  the  availability 
of  a  kanban  ticket  in  place  p*a„6an,  .  Then,  the  pallet  proceeds  to  the  machine,  in  place  pmi.  After  being 
worked  by  L„lt,  a  part  is  checked  for  quality  and  it  is  either  transported  back  to  pm.  by  tbacki  for  further 
rework,  or  moved  out  of  the  machine  by  touti.  The  numerical  values  of  the  parameters  for  the  model  are: 
w„M  =  1.2,  w„l2  =  1.4,  w„l3  =  1.3,  and  w,„4  =  1.1,  while,  for  each  sub-GSPN  i,  woki  =  0.7,  wre<<0|  =  0.3, 
and  Whack,  -  0.3.  All  transitions  have  single-server  semantics.  The  input  and  output  rates  for  the  entire 
kanban  network  are  set  by  assigning  wiui  =  1.0  and  wouti  =  0.9.  The  synchronizing  transition  tsyndli 
corresponds  to  the  merging  of  transitions  and  has  rate  0.4,  while  tsynch2  corresponds 

to  transitions  {tout2 ,  toul3 ,  t,Ui }  and  has  rate  0.5.  In  the  computations,  we  vary  the  number  N  =  N{  of 
tokens  initially  in  each  place  ])kunban, ■ 

We  consider  two  alternative  decompositions  of  the  model.  In  Case  1 ,  each  node  of  the  kanban  network 
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corresponds  to  a  sub-GSPN,  or  one  macroplace,  in  our  terminology.  Hence,  [P\  =  4  and  |A'#|  =  2.  In 
this  ca se,  the  fork-join  synchronization  creates  unreachable  markings  (any  marking  where  the  number 
of  tokens  in  Pkanban2  and  pkanban3  differ),  so  that  |Sr|  <  \S'\.  Instead,  in  Case  2,  we  merge  sub-GSPNs  2 
and  3  into  a  single  macroplace,  hence  \V\  =  3.  Now,  however,  the  constraint  on  the  number  of  tokens 
in  pkanban2  and  Pkanbanz  is  taken  explicitly  into  account  when  generating  the  local  state  space  for  the 
corresponding  sub-GSPN.  Thus,  the  overall  state  space  becomes  the  exact  cross-product  of  the  three 
local  state  spaces,  \St\  ==  |<S'| ,  and  the  computation  is  more  efficient,  since  there  is  no  need  to  distinguish 
the  reachable  states  from  the  unreachable  ones.  This  shows  how  an  intelligent  decomposition  of  the 
model  can  greatly  affect  the  computational  requirements  of  the  solution. 

Finally,  we  make  the  two  synchronizing  transitions  immediate  and  apply  Theorem  5.2  to  a  decom¬ 
position  into  either  four  or  three  sub-GSPNs,  resulting  in  Case  3  and  4,  respectively.  We  observe  that 
the  presence  of  synchronizing  transitions  leads  to  unreachable  vanishing  markings  in  this  model. 


N 

Case 

el 

e2 

64 

r 

1 

1,2 

0.907 

0.671 

0.671 

0.355 

0.132 

2 

1.810 

1.328 

1.328 

0.764 

0.248 

3 

2.722 

1.944 

1.944 

1.154 

0.332 

4 

3.646 

2.515 

2.515 

1.510 

0.393 

5 

4.588 

3.053 

3.053 

1.876 

0.439 

1 

3,4 

0.860 

0.810 

0.810 

0.534 

0.139 

2 

1.691 

1.671 

1.671 

1.268 

0.263 

3 

2.529 

2.545 

2.545 

2.037 

0.348 

4 

3.379 

3.425 

3.425 

2.807 

0.409 

5 

3.789 

4.454 

4.454 

3.914 

0.529 

Table  1:  Performance  results  as  a  function  of  N 

We  compute  the  expected  number  e:  of  tokens  in  places  pm. ,  pbackn  and  poult  for  each  sub-GSPN  of 
the  kanban  system.  For  a  given  i  =  1,...,4,  this  corresponds  to  a  reward  structure  where  the  reward 
rate  p( m)  of  a  marking  m  is  given  by  mm|.  +  m&adt-<  +  mouti)  and  the  reward  impulses  r  are  identically 
zero.  We  also  compute  the  throughput  r  of  the  system,  defined  as  the  expected  firing  rate  of  tsynchx  •  This 
corresponds  to  a  reward  structure  where  the  reward  rates  are  identically  zero  and  the  reward  impulses 
are  one  if  they  correspond  to  the  firing  of  tsynchix  zero  otherwise  (the  same  quantity  could  be  computed 
by  observing  tini ,  tsynch2 ,  or  touti  instead).  We  stress  that,  in  general,  reward  rates  and  impulses  could  be 
marking-dependent.  Table  1,  shows  the  resulting  numerical  values;  as  expected,  e2  =  e3.  It  is  apparent 
that  the  first  sub-GSPN  is  the  most  loaded  for  Case  1  and  2,  while  the  second  and  third  sub-GSPNs  are 
the  bottleneck  in  Case  3  and  4.  This  would  suggest  that,  the  synchronization  behavior  is  an  important 
performance  factor,  in  addition  to  the  actual  machining  rate  of  each  station. 

Table  2  shows  the  size  of  the  state  spaces  St ,  <Sy,  and  S'  as  a  function  of  N.  The  overall  number 
of  nonzero  elements  for  the  sparse  matrices  used  by  the  Kronecker  approach  (“local”)  can  be  compared 
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N 

Case 

|<ST| 

\Sv\ 

i«n 

local 

nonzero(Q) 

nonzero(P) 

“explore55 

“solve” 

1 

1 

160 

0 

256 

20 

616 

616 

0.017 

2 

4,600 

0 

10,000 

80 

28,120 

28,120 

0.967 

3 

58,400 

0 

160,000 

200 

446,400 

446,400 

16.267 

309.533 

4 

454,475 

0 

1,500,625 

400 

3,979,850 

3,979,850 

190.500 

4,721.217 

5 

2,546,432 

0 

9,834,496 

700 

24,460,016 

24,460,016 

1,759.000 

22,215.257 

1 

2 

160 

0 

160 

40 

616 

616 

not  req. 

0.017 

2 

4,600 

0 

4,600 

186 

28,120 

28,120 

not  req. 

1.517 

3 

58,400 

0 

58,400 

678 

446,400 

446,400 

not  req. 

50.483 

4 

454,475 

0 

454,475 

1878 

3,979,850 

3,979,850 

not  req. 

855.917 

5 

2,546,432 

0 

2,546,432 

4368 

24,460,016 

24,460,016 

not  req. 

6,054.783 

1 

3 

152 

8 

256 

20 

600 

608 

0.017 

0.167 

2 

3,816 

697 

10,000 

80 

23,832 

24,529 

0.517 

20.450 

3 

41,000 

13,656 

160,000 

200 

316,360 

330,016 

5.967 

84.600 

4 

268,475 

128,000 

1,500,625 

400 

2,343,050 

2,471,050 

75.783 

719.050 

5 

1,270,962 

769,480 

9,834,496 

700 

12,025,566 

12,795,046 

707.483 

4,111.681 

1 

4 

152 

8 

160 

40 

600 

608 

0.013 

0.150 

2 

3,816 

697 

4,600 

186 

23,832 

24,529 

0.333 

18.550 

3 

41,000 

13,656 

58,400 

678 

316,360 

330,016 

5.550 

76.250 

4 

268,475 

128,000 

454,475 

1878 

2,343,050 

2,471,050 

67.750 

647.100 

5 

1,270,962 

769,480 

2,546,432 

4368 

12,025,566 

12,795,046 

725.333 

3,562.531 

Table  2:  Computational  and  storage  requirements  as  a  function  of  N. 

to  that  of  conventional  solution  methods  that  explicitly  generate  Q.  For  comparison,  we  also  list  the 
number  of  nonzero  elements  in  P,  that  is  after  eliminating  the  local  vanishing  markings.  For  Case  1  and 
2,  this  makes  no  difference  (there  are  no  synchronizing  vanishing  markings),  while,  for  Case  3  and  4,  the 
memory  requirements  are  somewhat  larger.  It  is  clear  that  the  only  practical  limitation  of  our  approach 
is  the  memory  required  by  the  two  iteration  vectors,  and  7rlm+1l,  and  the  addressing  vector  $ j 
(of  sizes  |«St|),  or  by  7^,  and  ^  (of  sizes  \ST\  +  |<Sy|)?  plus  the  vector  h  (of  size  |<Sr|),  if  not 

computed  at  each  iteration.  It  is  also  apparent  how  the  approach  allows  the  solution  of  problems  one 
order  of  magnitude  larger  than  with  conventional  methods  in  an  acceptable  amount  of  time. 

The  computation  times  “explore55,  to  explore  the  reachability  set,  and  “solve55,  for  the  Jacobi  method 
with  relaxation  parameter  d  =  0.9,  are  given  in  seconds.  The  vector  h  is  stored  explicitly.  The 
convergence  criterion  is  set  to  ||7rtml  —  7rlm+1l||oo  <  10‘6  (or  H7M  —  7fm+1l||oo  <  10“6).  The  uniform 
distribution  was  used  as  the  initial  guess  for  7 (or  7I0]).  The  program  was  run  on  a  Sony  NWS-5000 
workstation  with  90  Mbyte  of  main  memory. 

We  stress  that  the  elimination  approach  (Theorem  4.1)  should  be  used  whenever  possible,  and  that 
preservation  of  the  non-local  vanishing  markings  (Theorem  5.2)  should  be  used  only  when  there  are 
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immediate  synchronizing  transitions.  Preservation  of  all  vanishing  markings  (Theorem  5.1)  is  probably 
never  appropriate,  since  memory  usage  is  the  paramount  consideration,  and  this  approach  increases  the 
already  critical  size  of  the  probability  vectors.  In  few  pathological  cases,  it  could  reduce  the  size  of  the 
“local”  data  structures,  but  these  are  not  critical. 

7  Conclusion 

We  rigorously  formalized  and  implemented  an  approach  based  on  Kronecker  algebra  for  the  solution  of 
the  CTMC  underlying  a  GSPN.  The  results  extend  previous  works  by  Donatelli,  Buchholz,  and  Kem¬ 
per  [4,  5,  10,  11,  12],  to  include  immediate  synchronizing  transitions,  quite  general  marking-dependent 
behavior,  and  a  reward  structure  allowing  reward  impulses  associated  with  immediate  transition  firings. 
The  restrictions  imposed  on  the  GSPN  are  minimal,  thus  the  approach  has  obvious  practical  appli¬ 
cations.  Furthermore,  the  structure  of  the  GSPN  itself  gives  strong  hints  on  its  decomposition.  For 
example,  if  an  “elimination-based”  solution  is  desired,  the  GSPN  must  be  decomposed  so  that  immedi¬ 
ate  transitions  are  local  to  a  sub-GSPN;  if  the  marking-dependency  of  a  transition  does  not  satisfy  our 
requirements,  all  the  places  responsible  for  this  behavior  should  also  be  merged  in  the  same  sub-GSPN. 

Memory  requirements  are  still  the  main  limitation  to  the  solution,  but  these  have  now  been  reduced 
from  the  size  of  the  transition  rate  matrix  to  that  of  the  steady-state  probability  vector,  for  a  very  general 
class  of  GSPNs.  Even  for  highly  sparse  matrices,  this  corresponds  to  the  ability  to  solve  problems  whose 
state  space  is  one  order  of  magnitude  larger  than  with  a  traditional  solution  approach. 

To  further  increase  the  size  of  models  that  can  be  solved,  we  can  use  the  distributed  state- space 
generation  algorithm  described  in  [7],  which  allows  one  to  partition  the  memory  and  execution  require¬ 
ments  to  generate  the  state  space  over  a  set  of  workstations,  and  which  exhibits  excellent  speedups  for 
large  problems.  The  Jacobi  method  we  employed  can  also  be  parallelized  using  a  set  of  workstations,  so 
that  the  entire  solution  process  is  performed  in  a  distributed  fashion  and  uses  the  available  memory.  In 
particular,  the  “local”  matrices  occupy  a  negligible  amount  of  memory  and  can  be  duplicated  on  each 
processor,  while  only  the  probability  vector  needs  to  be  distributed. 

Since  our  results  are  complementary  to  those  in  [7],  we  can  reasonably  hope  for  a  two-orders  of  mag¬ 
nitude  increase  in  the  size  of  manageable  models,  assuming  the  availability  of  a  network  of  workstations. 
As  a  target  for  the  near  future,  we  will  attempt  to  solve  CTMCs  with  10s  states  and  a  transition  rate 
matrix  with  109  non-zeros  in  a  matter  of  hours  using  a  dozen  workstations  with  128Mbytes  of  memory 
each. 

Eventually,  even  CTMCs  of  this  size  might  not  be  enough  to  satisfy  the  needs  of  a  serious  modeler. 
Our  results,  combined  with  the  distributed  state-space  generation  algorithm,  are  directly  applicable  to 
the  family  of  approximations  suggested  in  [17].  We  can  simply  consider  each  sub-GSPN  separately  for 
the  coarsest  model,  two  adjacent  sub-GSPNs  for  the  next  more  accurate  model,  then  three  adjacent 
sub-GSPNs  and  so  on.  The  comparison  of  the  results  from  successive  approximate  models  can  suggest 
an  estimate  of  the  quality  of  the  results.  This  allows  trading  off  lower  approximation  errors  against 
higher  computational  effort.  We  are  currently  investigating  the  requirements  to  assure  an  accurate 
solution  (e.g.,  approximation  errors  below  5%)  on  a  single  workstation  for  a  large  class  of  models,  where 
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the  CTMCs  may  have  more  than  1030  states. 

A  prototype  version  of  the  program  used  to  compute  the  results  we  presented  is  available  and  can 
be  obtained  by  contacting  the  second  author.  The  input  to  the  program  has  a  SPNP-style  syntax  [9]. 
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